Plotting invariant sets for periodic type wind-tree models¶
from surface_dynamics import *
import sympy as sp
RF = RealField(256)
Find periodic parameters for genus 2 surface $Y$¶
We find a periodic type IET corresponding to the flow on $Y$ by manually taking a loop in the Rauzy graph. Then we need to check whether this in fact corresponds to some parameters $(a,b,t)$.
perm_Y = iet.Permutation('a b c d e','e a c b d')
d = perm_Y.rauzy_diagram()
d.graph().show(color_by_label=true,figsize=30,layout='spring',iterations=200)
perm_Y = iet.Permutation('a b c d e','e a c b d')
path = ['t', 't', 'b', 't', 'b', 't', 't','b','t','b','t','b','b','t']
A = iet.FlipSequence(perm_Y,path).matrix()
A = matrix(SR,A)
eig = A.eigenvectors_right()
eig.sort(key=lambda eigenspace: eigenspace[0])
maxevec = eig[-1][1][0]
maxevec = maxevec/sum(maxevec)
param_a = RF(maxevec[1])
param_t = RF(1/(maxevec[4]+1-param_a))
param_b = RF(param_t*(-maxevec[2]+(1-param_a)))
Construct IET for the flow on $X$ with parameters found above¶
#Define the linear flow on the genus 5 fourfold cover surface
#Denote the four copies of the genus 2 surface by (0,0), (1,0), (0,1), (0,0)
#Start from point (x,y) in a given copy and go forwards until intersection with boundary
a = param_a
b = param_b
t = param_t
def forward(x,y,copy):
#Find the first intersections with the four possible sides, variable "side" is the
#horizontal distance until intersection
#We also record the number of intersections with outer edges, this will be useful for
#finding the cocycle determining the cover
intersections = {}
for c in [(0,0),(1,0),(0,1),(1,1)]:
intersections[('h',c)] = 0
intersections[('v',c)] = 0
#inner vertical
if x < (1-a)/2:
#check really hit inner edge rather than its extension
innerv = (1-a)/2 - x
if y+innerv*t >= (1-b)/2 and y+innerv*t <= (1+b)/2:
pass
else:
innerv = Infinity
else:
innerv = Infinity
#inner horizontal
if y < (1-b)/2:
#check really hit inner edge rather than its extension
innerh = ((1-b)/2 - y)/t
if x+innerh >= (1-a)/2 and x+innerh <= (1+a)/2:
pass
else:
innerh = Infinity
else:
innerh = Infinity
#outer vertical
outerv = 1-x
#outer horizontal
outerh = (1-y)/t
#find minimum to see which intersection happens
if min(innerv,innerh,outerv,outerh) == innerv:
hitx = (1-a)/2
hity = y+innerv*t
newx = (1+a)/2 #jump to glued edge
newy = y+innerv*t
newcopy = ((copy[0]+1) % 2, copy[1])
if min(innerv,innerh,outerv,outerh) == innerh:
hitx = x+innerh
hity = (1-b)/2
newx = x+innerh
newy = (1+b)/2
newcopy = (copy[0], (copy[1]+1) % 2)
if min(innerv,innerh,outerv,outerh) == outerv:
hitx = 1
hity = y+outerv*t
newx = 0
newy = y+outerv*t
newcopy = copy
intersections[('v',copy)] += 1
if min(innerv,innerh,outerv,outerh) == outerh:
hitx = x+outerh
hity = 1
newx = x+outerh
newy = 0
newcopy = copy
intersections[('h',copy)] += 1
#return 2 points: before and after gluing
return((hitx,hity,copy),(newx,newy,newcopy),intersections)
def backward(x,y,copy):
#Flowing backwards
#inner vertical
if x > (1+a)/2:
#check really hit inner edge rather than its extension
innerv = x - (1+a)/2
if y-innerv*t >= (1-b)/2 and y-innerv*t <= (1+b)/2:
pass
else:
innerv = Infinity
else:
innerv = Infinity
#inner horizontal
if y > (1+b)/2:
#check really hit inner edge rather than its extension
innerh = (y - (1+b)/2)/t
if x-innerh >= (1-a)/2 and x-innerh <= (1+a)/2:
pass
else:
innerh = Infinity
else:
innerh = Infinity
#outer vertical
outerv = x
#outer horizontal
outerh = (y)/t
#find minimum to see which intersection happens
if min(innerv,innerh,outerv,outerh) == innerv:
hitx = (1+a)/2
hity = y-innerv*t
newx = (1-a)/2 #jump to glued edge
newy = y-innerv*t
newcopy = ((copy[0]+1) % 2, copy[1])
if min(innerv,innerh,outerv,outerh) == innerh:
hitx = x-innerh
hity = (1+b)/2
newx = x-innerh
newy = (1-b)/2
newcopy = (copy[0], (copy[1]+1) % 2)
if min(innerv,innerh,outerv,outerh) == outerv:
hitx = 0
hity = y-outerv*t
newx = 1
newy = y-outerv*t
newcopy = copy
if min(innerv,innerh,outerv,outerh) == outerh:
hitx = x-outerh
hity = 0
newx = x-outerh
newy = 1
newcopy = copy
#return 2 points: before and after gluing
return((hitx,hity,copy),(newx,newy,newcopy))
def section_fwd(x,y,copy):
#flow forward until Poincare section, recording the number of interseections with edges
intersections = {}
for c in [(0,0),(1,0),(0,1),(1,1)]:
intersections[('h',c)] = 0
intersections[('v',c)] = 0
((x,y,copy),step_intersections) = forward(x,y,copy)[1:]
for curve in intersections.keys():
intersections[curve] += step_intersections[curve]
while True:
if y == 0 and copy == (0,0):
return (x,intersections)
else:
((x,y,copy),step_intersections) = forward(x,y,copy)[1:]
for curve in intersections.keys():
intersections[curve] += step_intersections[curve]
def section_back(x,y,copy):
#flow backward until Poincare section
(x,y,copy) = backward(x,y,copy)[1]
while True:
if y == 1 and copy == (0,0):
return x
else:
(x,y,copy) = backward(x,y,copy)[1]
Find the IET
#Lengths
#Find the forwards discontinuities of the IET by flowing backwards from the corner
discont = []
vertices = [((1-a)/2,(1-b)/2),((1+a)/2,(1-b)/2),((1-a)/2,(1+b)/2)]
for i in range(2):
for j in range(2):
copy = (i,j)
for v in vertices:
(x,y) = v
discont.append(section_back(x,y,copy))
discont.append(section_back(1,1,(0,0))) #end of section
discont.sort()
discont = [0] + discont + [1]
lengths0 = [discont[i+1] - discont[i] for i in range (len(discont)-1)]
#Permutation
#for each discontinuity (including 0 but not 1) find its right image
l = min(lengths0)
discont_right_images = [section_fwd(x+l/2,0,(0,0))[0]-(l/2) for x in discont[:-1]]
#Now the permutation is just the order of these images
letters = ['a','b','c','d','e','f','g','h','i','j','k','l','m','n']
unordered = list(zip(letters,discont_right_images))
ordered = sorted(unordered,key = lambda pt: pt[1])
permuted = [letter[0] for letter in ordered]
perm0 = iet.Permutation(letters,permuted)
T0 = iet.IntervalExchangeTransformation(perm0,lengths0)
Find the cocycles $\gamma_h, \gamma_v$ corresponding to the cover
gamma_h = []
gamma_v = []
sigma_h = []
sigma_v = []
for x in discont[:-1]:
intersections = section_fwd(x+l/2,0,(0,0))[1]
gamma_h.append(intersections[('v',(0,0))]-intersections[('v',(1,0))]+intersections[('v',(0,1))]-intersections[('v',(1,1))])
gamma_v.append(intersections[('h',(0,0))]+intersections[('h',(1,0))]-intersections[('h',(0,1))]-intersections[('h',(1,1))])
sigma_v.append(intersections[('v',(0,0))]+intersections[('v',(1,0))]-intersections[('v',(0,1))]-intersections[('v',(1,1))])
sigma_h.append(intersections[('h',(0,0))]-intersections[('h',(1,0))]+intersections[('h',(0,1))]-intersections[('h',(1,1))])
gamma_h = vector(gamma_h)
gamma_v = vector(gamma_v)
sigma_h = vector(sigma_h)
sigma_v = vector(sigma_v)
Find matrix corresponding to the period of Rauzy-Veech induction for $T$¶
T0 = iet.IntervalExchangeTransformation(perm0,lengths0)
T = T0
path = []
for i in range(1,10000):
(T,move) = T.rauzy_move(data=True)[:2]
path.append(move)
if T.permutation() == perm0:
print("Found loop of length "+str(i))
break
A = matrix(iet.FlipSequence(perm0,path).matrix())
Found loop of length 434
In fact taking half the loop we obtain the same permutation up to relabeling
T0 = iet.IntervalExchangeTransformation(perm0,lengths0)
T = T0
path0 = []
for i in range(217):
(T,move) = T.rauzy_move(data=True)[:2]
path0.append(move)
print(T.permutation(),"\n",perm0)
a k l g h i d e f m b c j n n a c b m l e d i h g f k j a b c d e f g h i j k l m n n a l k j c h g f e d i b m
#Check that indeed just have a relabeling
B = matrix(SR,iet.FlipSequence(perm0,path0,relabelling=[0,10,11,6,7,8,3,4,5,12,1,2,9,13]).matrix())
B^2 == A
True
#Now that we have the matrix e can compute the lengths exactly
Bsp = sp.Matrix(B)
eig = Bsp.eigenvects()
eig.sort(key = lambda espace: espace[0])
lens_evalue = eig[-1][0]
precise_lengths0 = vector(eig[-1][2][0])
precise_lengths0 /= sum(precise_lengths0)
Find stable vectors¶
Now that we have the matrix $B$, we can find the stable vectors in the blocks $E^{-+},E^{+-}$ by taking the eignenvectors of $B$. To simplify, we only consider the action of $B$ on the individual blocks, which also allows us to write the stable vectors in the basis $\{\gamma,\sigma\}$ for the corresponding block.
#Find B^T * gamma_h in terms of gamma_h, sigma_h
#orthogonal vector to gamma_h
gamma_h_orth = sigma_h - (sigma_h*gamma_h)/(gamma_h*gamma_h) *gamma_h
#project onto orthogonal vector to find the sigma_h components, then rest is gamma_h component
image_gamma_h = B.transpose()*gamma_h
c = (image_gamma_h*gamma_h_orth)/(sigma_h*gamma_h_orth)
a = (image_gamma_h-c*sigma_h)*gamma_h/(gamma_h*gamma_h)
image_sigma_h = B.transpose()*sigma_h
d = (image_sigma_h*gamma_h_orth)/(sigma_h*gamma_h_orth)
b = (image_sigma_h-d*sigma_h)*gamma_h/(gamma_h*gamma_h)
#matrix for action of B on the horizontal block
h_matrix = matrix(SR,[[a,b],[c,d]])
#check that the below are equal to 0
# image_gamma_h - h_matrix[0][0]*gamma_h - h_matrix[1][0]*sigma_h == 0
# image_sigma_h - h_matrix[0][1]*gamma_h - h_matrix[1][1]*sigma_h == 0
h_eigenvecs = h_matrix.eigenvectors_right()
h_eigenvecs.sort(key= lambda evec: evec[0])
h_const = (h_eigenvecs[0][1][0][1]/h_eigenvecs[0][1][0][0]) #the sigma entry of the stable eigenvector
h_stable = gamma_h + h_const*sigma_h
#lambda, the scaling factor
lam = h_eigenvecs[0][0]
print(lam)
#check that really stable eigenvector:
[entry.full_simplify() for entry in (B.transpose()*h_stable - lam*h_stable)]
-20*sqrt(6) + 49
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
#Find B^T * gamma_v in terms of gamma_v, sigma_v
#orthogonal vector to gamma_v
gamma_v_orth = sigma_v - (sigma_v*gamma_v)/(gamma_v*gamma_v) *gamma_v
#project onto orthogonal vector to find the sigma_v components, then rest is gamma_v component
image_gamma_v = B.transpose()*gamma_v
c = (image_gamma_v*gamma_v_orth)/(sigma_v*gamma_v_orth)
a = (image_gamma_v-c*sigma_v)*gamma_v/(gamma_v*gamma_v)
image_sigma_v = B.transpose()*sigma_v
d = (image_sigma_v*gamma_v_orth)/(sigma_v*gamma_v_orth)
b = (image_sigma_v-d*sigma_v)*gamma_v/(gamma_v*gamma_v)
#matrix for action of B on the horizontal block
v_matrix = matrix(SR,[[a,b],[c,d]])
#check that the below are equal to 0
# image_gamma_v - v_matrix[0][0]*gamma_v - v_matrix[1][0]*sigma_v == 0
# image_sigma_v - v_matrix[0][1]*gamma_v - v_matrix[1][1]*sigma_v == 0
v_eigenvecs = v_matrix.eigenvectors_right()
v_eigenvecs.sort(key= lambda evec: evec[0])
v_const = v_eigenvecs[0][1][0][1]/v_eigenvecs[0][1][0][0] #the sigma entry of the stable eigenvector
v_stable = gamma_v + v_const*sigma_v
#check that really stable eigenvector:
[entry.full_simplify() for entry in (B.transpose()*v_stable - lam*v_stable)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Transfer functions via Birkhoff sums¶
Here we plot the transfer functions to see what they look like. For each of the two, we plot the transfer function $h$ and its translate by $c$
#define piecewise constant function corresponding to stable vector
def psi_h(x):
for i in range(14):
if x < discont[i+1]:
return h_stable[i]
def psi_v(x):
for i in range(14):
if x < discont[i+1]:
return v_stable[i]
N = 10^4
#values of the transfer functions
x0 = 0
current_h = 0
current_v = 0
x = x0
values_h = [(x,current_h)]
values_v = [(x,current_v)]
for i in range(N):
current_h += psi_h(x)
current_v += psi_v(x)
x = T0(x)
values_h.append((x,current_h))
values_v.append((x,current_v))
h_plt = points(values_h,figsize=100,aspect_ratio=0.01)+points([(val[0],val[1]+h_const) for val in values_h],color='red')
h_plt
v_plt = points(values_v,figsize=100,aspect_ratio=0.01)+points([(val[0],val[1]+v_const) for val in values_v],color='red')
v_plt
Transfer functions at endpoints of intervals¶
Here we find the vectors $\tau$ for each of the two transfer functions, defined as $(\tau)_i = h (x_i) - h (x_{i-1})$, where $x_i$ are the discontinuities of $T$.
We know that the vector $\tau$ is an eigenvector of B with eigenvalue $\lambda^{-1}$. (See Lemma A.1)
The eigenspace with eignavlue $\lambda^{-1}$ is 2-dimensional, so we need to identify our vectors inside this eigenspace.
Consider the permutation and observe that $\pi_b(1) = 14, \pi_b(14) = 13$. We know that if $\tau_h$ is the eigenvector for the horizontal stable cocycle, then:
$\sum_{i=1}^{13} (\tau_h)_i + h_{14} = h(0) = 0$,
$\sum_{i=1}^{13} (\tau_h)_i + h_{13} = h(1) = \sum_{i=1}^{14} (\tau_h)_i \implies (\tau_h)_{14} = h_{13}$.
Writing in a basis $\tau_1, \tau_2$ of the eigenspace, let $s_1 = \sum_{i=1}^{13} (\tau_1)_i, s_2 = \sum_{i=1}^{13} (\tau_2)_i$.
Then if $\tau_h = a \tau_1 + b \tau_2$, we have: $$ \begin{pmatrix} s_1 & s_2 \\ (\tau_1)_{14} & (\tau_2)_{14} \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} -h_{14} \\ h_{13} \end{pmatrix}. $$
Bsp = sp.Matrix(B)
eig = Bsp.eigenvects()
for espace in eig:
if SR(espace[0])*lam == 1 :
[tau1,tau2] = espace[2]
break
s_1 = sum(tau1[:13])
s_2 = sum(tau2[:13])
M = matrix(SR,[[s_1,s_2],[tau1[13],tau2[13]]])
(a,b) = M.inverse()*vector([-h_stable[13],h_stable[12]])
tau_h = a*tau1+b*tau2
(c,d) = M.inverse()*vector([-v_stable[13],v_stable[12]])
tau_v = c*tau1+d*tau2
tau_h = vector(SR,tau_h)
tau_v = vector(SR,tau_v)
Compute max/min values of $h$ for level 0 partition¶
We use the formula of Proposition 5.3 to find upper and lower bounds of $h$ on each $I_i$.
g = iet.FlipSequence(perm0,path0,relabelling=[0,10,11,6,7,8,3,4,5,12,1,2,9,13])
orbit_sub = g.orbit_substitution
orbit_sub_list = orbit_sub(' ')
interval_sub = g.interval_substitution
interval_sub_list = interval_sub(' ')
#values of the function f_h
#this we do with reals with high precision for speed, symbolic computation is much slower
f_h_values = []
#list, with entry t being {s:[values for edges s->t]}
for t in range(14):
vals = {i:[] for i in range(14)}
current = 0
for s in orbit_sub_list[t]:
vals[s].append(RF(current))
current += h_stable[s]
f_h_values.append(vals)
#max and min of f depending on start and terminus
h_maxmatrix = matrix([[max(f_h_values[j][i]) for j in range(14)] for i in range(14)])
h_minmatrix = matrix([[min(f_h_values[j][i]) for j in range(14)] for i in range(14)])
total_max = max([h_maxmatrix[i][j] for j in range(14) for i in range(14)])
total_min = min([h_minmatrix[i][j] for j in range(14) for i in range(14)])
print("Extrema for f: "+str((total_min,total_max)))
Extrema for f: (-1.974291367869850351808974451293155894823188058130298334761510844497815838545, 0.9844965122062883878632929571753280555042384449968957661076594994786082893959)
We can see that the max absolute value is less than 2, and since $\lambda \approx 0.01$, max value times $\frac{\lambda^4}{1-\lambda}$ is less that $3\times10^{-8}$.
So we compute the max and min for four steps and add the error afterwards.
def maxplusproduct(A,B):
#given matrices A and B return AB in max-plus arithmetic
n = len(list(A))
rows = []
for i in range(n):
row = []
for j in range(n):
row.append(max([A[i][k]+B[k][j] for k in range(n)]))
rows.append(row)
return matrix(rows)
def minplusproduct(A,B):
#given matrices A and B return AB in min-plus arithmetic
n = len(list(A))
rows = []
for i in range(n):
row = []
for j in range(n):
row.append(min([A[i][k]+B[k][j] for k in range(n)]))
rows.append(row)
return matrix(rows)
h_max2 = maxplusproduct(h_maxmatrix,lam*h_maxmatrix)
h_max3 = maxplusproduct(h_max2,lam^2*h_maxmatrix)
h_max4 = maxplusproduct(h_max3,lam^3*h_maxmatrix)
h_max_vec = [max(h_max4[i]) for i in range(14)]
h_min2 = minplusproduct(h_minmatrix,lam*h_minmatrix)
h_min3 = minplusproduct(h_min2,lam^2*h_minmatrix)
h_min4 = minplusproduct(h_min3,lam^3*h_minmatrix)
h_min_vec = [max(h_min4[i]) for i in range(14)]
#values of the function f_v
f_v_values = []
#list, with entry t being {s:[values for edges s->t]}
for t in range(14):
vals = {i:[] for i in range(14)}
current = 0
for s in orbit_sub_list[t]:
vals[s].append(RF(current))
current += v_stable[s]
f_v_values.append(vals)
v_maxmatrix = matrix([[max(f_v_values[j][i]) for j in range(14)] for i in range(14)])
v_minmatrix = matrix([[min(f_v_values[j][i]) for j in range(14)] for i in range(14)])
total_max = max([v_maxmatrix[i][j] for j in range(14) for i in range(14)])
total_min = min([v_minmatrix[i][j] for j in range(14) for i in range(14)])
print("Extrema for f: "+str((total_min,total_max)))
Extrema for f: (-0.05643925929693713551019126279916619925695797905930175009814259165278802838182, 0.7709372431646665952379706105474491594426470186321003481763863077354391648243)
v_max2 = maxplusproduct(v_maxmatrix,lam*v_maxmatrix)
v_max3 = maxplusproduct(v_max2,lam^2*v_maxmatrix)
v_max4 = maxplusproduct(v_max3,lam^3*v_maxmatrix)
v_max_vec = [max(v_max4[i]) for i in range(14)]
v_min2 = minplusproduct(v_minmatrix,lam*v_minmatrix)
v_min3 = minplusproduct(v_min2,lam^2*v_minmatrix)
v_min4 = minplusproduct(v_min3,lam^3*v_minmatrix)
v_min_vec = [max(v_min4[i]) for i in range(14)]
Get bounds on partition¶
Compute the relative min/max of $h$ on interval $I_i$ compared to the left endpoint: the $i$th entry of h_relmax is $\max_{y\in I_i} h(y) - h(x_i)$. (Up to the error of size $3\times 10^{-8}$)
h_relmax = [float(h_max_vec[i]-sum(tau_h[:i])) for i in range(14)]
h_relmin = [float(h_min_vec[i]-sum(tau_h[:i])) for i in range(14)]
v_relmax = [float(v_max_vec[i]-sum(tau_v[:i])) for i in range(14)]
v_relmin = [float(v_min_vec[i]-sum(tau_v[:i])) for i in range(14)]
Now we compute the min/max bounds on all intervals of the level one partition $\mathcal{P}^{(1)}$. This takes some time! (ca 20 minutes on my computer)
err = 3*10^(-8)
x = 0
current_h = 0
current_v = 0
h_bounds = []
v_bounds = []
for i in range(14):
print("Computing on interval I_"+str(i))
for interval in interval_sub_list[i]:
h_bounds.append((float(x),float(current_h+lam*h_relmin[interval]-err),float(current_h+lam*h_relmax[interval]+err)))
v_bounds.append((float(x),float(current_v+lam*v_relmin[interval]-err),float(current_v+lam*v_relmax[interval]+err)))
x += 1/lens_evalue*precise_lengths0[interval]
current_h += lam*tau_h[interval]
current_v += lam*tau_v[interval]
Computing on interval I_0 Computing on interval I_1 Computing on interval I_2 Computing on interval I_3 Computing on interval I_4 Computing on interval I_5 Computing on interval I_6 Computing on interval I_7 Computing on interval I_8 Computing on interval I_9 Computing on interval I_10 Computing on interval I_11 Computing on interval I_12 Computing on interval I_13
Now we split the intervals into groups and cmpute min/max over each, for faster comparisons later. We do this on two levels.
int_n = len(h_bounds)
bin_n = floor(int_n^(1/3))
ints_per_bin_1 = int_n//(bin_n^2)
#level 1 partition
bin_limits_1 = [(i*ints_per_bin_1,(i+1)*ints_per_bin_1) for i in range(bin_n^2-1)]+[((bin_n^2-1)*ints_per_bin_1,int_n-1)]
#level 2 partition
ints_per_bin_2 = len(bin_limits_1)//bin_n
bin_limits_2 = [(i*ints_per_bin_2,(i+1)*ints_per_bin_2) for i in range(bin_n-1)]+[((bin_n-1)*ints_per_bin_2,len(bin_limits_1)-1)]
bin_1_h_bounds = []
bin_1_v_bounds = []
for i in range(bin_n^2):
(start,end) = bin_limits_1[i]
min_h = min([bound[1] for bound in h_bounds[start:end]])
max_h = max([bound[2] for bound in h_bounds[start:end]])
min_v = min([bound[1] for bound in v_bounds[start:end]])
max_v = max([bound[2] for bound in v_bounds[start:end]])
bin_1_h_bounds.append((min_h,max_h))
bin_1_v_bounds.append((min_v,max_v))
bin_2_h_bounds = []
bin_2_v_bounds = []
for i in range(bin_n):
(start,end) = bin_limits_2[i]
min_h = min([bound[0] for bound in bin_1_h_bounds[start:end]])
max_h = max([bound[1] for bound in bin_1_h_bounds[start:end]])
min_v = min([bound[0] for bound in bin_1_v_bounds[start:end]])
max_v = max([bound[1] for bound in bin_1_v_bounds[start:end]])
bin_2_h_bounds.append((min_h,max_h))
bin_2_v_bounds.append((min_v,max_v))
Finally we can define the function that we will use to find the level sets of the transfer functions
def check(a_h,a_v):
#return points x for which both h_x(x)-a_h lies in Z*c_h and similar for v
#first check the bins, then if in bin for both, check small intervals
validpts = []
validbins_1 = []
validbins_2 = []
for i in range(len(bin_2_h_bounds)):
(h_min,h_max) = bin_2_h_bounds[i]
test1 = float((h_min-a_h)/h_const)
test2 = float((h_max-a_h)/h_const)
if floor(test2) != floor(test1):
(v_min,v_max) = bin_2_v_bounds[i]
test1 = float((v_min-a_v)/v_const)
test2 = float((v_max-a_v)/v_const)
if floor(test2) != floor(test1):
validbins_2.append(i)
#check the small bins within the large bin
for i in validbins_2:
checkrange = bin_limits_2[i]
for j in range(checkrange[0],checkrange[1]+1):
(h_min,h_max) = bin_1_h_bounds[j]
test1 = float((h_min-a_h)/h_const)
test2 = float((h_max-a_h)/h_const)
if floor(test2) != floor(test1):
(v_min,v_max) = bin_1_v_bounds[j]
test1 = float((v_min-a_v)/v_const)
test2 = float((v_max-a_v)/v_const)
if floor(test2) != floor(test1):
validbins_1.append(j)
#check the intervals within the small bins
for i in validbins_1:
checkrange = bin_limits_1[i]
for j in range(checkrange[0],checkrange[1]+1):
(x,h_min,h_max) = h_bounds[j]
test1 = float((h_min-a_h)/h_const)
test2 = float((h_max-a_h)/h_const)
if floor(test2) != floor(test1):
(x,v_min,v_max) = v_bounds[j]
test1 = float((v_min-a_v)/v_const)
test2 = float((v_max-a_v)/v_const)
if floor(test2) != floor(test1):
validpts.append(x)
return validpts
Plot level sets¶
First we just make a plot of the obstacles in the windtree, and define the billiard map which we will use to plot segments of the trajectory between the hits of the Poincaré section
windtree_plt = plot([],figsize=200)
#make [-N,N]x[-N,N] grid
N = 90
a = param_a
b = param_b
for i in range(-N,N+1):
for j in range(-N,N+1):
windtree_plt+=polygon([(i-a/2,j-b/2),(i+a/2,j-b/2),(i+a/2,j+b/2),(i-a/2,j+b/2)],color='black',fill=True,alpha=1,axes=False)
def bil_map(pt,direction):
#image of billiard map: next reflection, next direction
#point given in coordinates, direction is of ray from horizontal
(x0,y0) = pt
#compute first intersection with each edge (compute the x-shift) - can be positive or negative
leftint = ceil(x0+a/2)-a/2-x0
botint = (ceil(y0+b/2)-b/2-y0)/tan(direction)
rightint = floor(x0-a/2)+a/2-x0
topint = (floor(y0-b/2)+b/2-y0)/tan(direction)
#find all intersections
#left
lshift = 0
n = 0
while lshift == 0:
shift = leftint+n
y = y0+tan(direction)*shift
if abs(round(y)-y)<b/2:
lshift = shift
break
n +=1
#bottom
bshift = 0
n = 0
while bshift == 0:
shift = botint+n/tan(direction)
x = x0+shift
if abs(round(x)-x)<a/2:
bshift = shift
break
n +=1
#right
rshift = 0
n = 0
while rshift == 0:
shift = rightint+n
y = y0+tan(direction)*shift
if abs(round(y)-y)<b/2:
rshift = shift
break
n -= 1
#top
tshift = 0
n = 0
while tshift == 0:
shift = topint+n/tan(direction)
x = x0+shift
if abs(round(x)-x)<a/2:
tshift = shift
break
n -= 1
if float(direction/(2*pi)) % 1 < 1/4:
#possible edges: left, bottom
if lshift < bshift:
shift = lshift
newdir = pi-direction
elif lshift > bshift:
shift = bshift
newdir = 2*pi-direction
elif float(direction/(2*pi)) % 1 < 2/4:
#possible edges: right, bottom
if rshift > bshift:
shift = rshift
newdir = pi-direction
elif rshift < bshift:
shift = bshift
newdir = 2*pi-direction
elif float(direction/(2*pi)) % 1 < 3/4:
#possible edges: right,top
if rshift > tshift:
shift = rshift
newdir = pi-direction
elif rshift < tshift:
shift = tshift
newdir = 2*pi-direction
elif float(direction/(2*pi)) % 1 < 4/4:
#possible edges: left, top
if lshift < tshift:
shift = lshift
newdir = pi-direction
elif lshift > tshift:
shift = tshift
newdir = 2*pi-direction
x = x0+shift
y = y0+tan(direction)*shift
return ((x,y),newdir)
def first_return(pt,direction):
#first hit of trajectory from pt in given direction to the section deifning IET
dir0 = direction
found = False
pts = [pt]
while not found:
(newpt,newdir) = bil_map(pt,direction)
if abs(float(direction/(2*pi))%1-float(dir0/(2*pi))%1) < 0.1:
#check if crossed horizontal line of form y = n+0.5
if floor(pt[1]+0.5) != floor(newpt[1]+0.5):
found = True
pts.append(newpt)
(pt,direction) = (newpt,newdir)
return pts
#Fix values of h to plot the level set
h0val = 0.2
v0val = -0.2
N=120
#Will look for the level set on the square [-N,N]x[-N,N]
#It's better to take this N larger than the N above, as around the edges we will miss parts of the
#trajectory that don't hit the Poincaré section inside the box
valid_pts = {}
for a_h in range(-N, N+1):
print("Computing column "+str(a_h))
for a_v in range(-N, N+1):
square = (a_h,a_v)
newvalid = []
for pt in check(a_h+h0val,a_v+v0val):
newvalid.append(pt)
if newvalid != []:
valid_pts[square] = newvalid
#To get the coordinates for plotting we need to shift by (-0.5,-0.5)
plot_pts = [(float(-0.5+x+square[0]),float(-0.5+square[1])) for square in valid_pts.keys() for x in valid_pts[square]]
theta = atan(param_t)
print("Plotting")
plt_copy = windtree_plt
for pt in plot_pts:
#plot first return to the section
traj = first_return(pt,float(theta))
for i in range(len(traj)-1):
plt_copy += line([traj[i],traj[i+1]],color='mediumblue',thickness=5)
plt_copy
Computing column -120 Computing column -119 Computing column -118 Computing column -117 Computing column -116 Computing column -115 Computing column -114 Computing column -113 Computing column -112 Computing column -111 Computing column -110 Computing column -109 Computing column -108 Computing column -107 Computing column -106 Computing column -105 Computing column -104 Computing column -103 Computing column -102 Computing column -101 Computing column -100 Computing column -99 Computing column -98 Computing column -97 Computing column -96 Computing column -95 Computing column -94 Computing column -93 Computing column -92 Computing column -91 Computing column -90 Computing column -89 Computing column -88 Computing column -87 Computing column -86 Computing column -85 Computing column -84 Computing column -83 Computing column -82 Computing column -81 Computing column -80 Computing column -79 Computing column -78 Computing column -77 Computing column -76 Computing column -75 Computing column -74 Computing column -73 Computing column -72 Computing column -71 Computing column -70 Computing column -69 Computing column -68 Computing column -67 Computing column -66 Computing column -65 Computing column -64 Computing column -63 Computing column -62 Computing column -61 Computing column -60 Computing column -59 Computing column -58 Computing column -57 Computing column -56 Computing column -55 Computing column -54 Computing column -53 Computing column -52 Computing column -51 Computing column -50 Computing column -49 Computing column -48 Computing column -47 Computing column -46 Computing column -45 Computing column -44 Computing column -43 Computing column -42 Computing column -41 Computing column -40 Computing column -39 Computing column -38 Computing column -37 Computing column -36 Computing column -35 Computing column -34 Computing column -33 Computing column -32 Computing column -31 Computing column -30 Computing column -29 Computing column -28 Computing column -27 Computing column -26 Computing column -25 Computing column -24 Computing column -23 Computing column -22 Computing column -21 Computing column -20 Computing column -19 Computing column -18 Computing column -17 Computing column -16 Computing column -15 Computing column -14 Computing column -13 Computing column -12 Computing column -11 Computing column -10 Computing column -9 Computing column -8 Computing column -7 Computing column -6 Computing column -5 Computing column -4 Computing column -3 Computing column -2 Computing column -1 Computing column 0 Computing column 1 Computing column 2 Computing column 3 Computing column 4 Computing column 5 Computing column 6 Computing column 7 Computing column 8 Computing column 9 Computing column 10 Computing column 11 Computing column 12 Computing column 13 Computing column 14 Computing column 15 Computing column 16 Computing column 17 Computing column 18 Computing column 19 Computing column 20 Computing column 21 Computing column 22 Computing column 23 Computing column 24 Computing column 25 Computing column 26 Computing column 27 Computing column 28 Computing column 29 Computing column 30 Computing column 31 Computing column 32 Computing column 33 Computing column 34 Computing column 35 Computing column 36 Computing column 37 Computing column 38 Computing column 39 Computing column 40 Computing column 41 Computing column 42 Computing column 43 Computing column 44 Computing column 45 Computing column 46 Computing column 47 Computing column 48 Computing column 49 Computing column 50 Computing column 51 Computing column 52 Computing column 53 Computing column 54 Computing column 55 Computing column 56 Computing column 57 Computing column 58 Computing column 59 Computing column 60 Computing column 61 Computing column 62 Computing column 63 Computing column 64 Computing column 65 Computing column 66 Computing column 67 Computing column 68 Computing column 69 Computing column 70 Computing column 71 Computing column 72 Computing column 73 Computing column 74 Computing column 75 Computing column 76 Computing column 77 Computing column 78 Computing column 79 Computing column 80 Computing column 81 Computing column 82 Computing column 83 Computing column 84 Computing column 85 Computing column 86 Computing column 87 Computing column 88 Computing column 89 Computing column 90 Computing column 91 Computing column 92 Computing column 93 Computing column 94 Computing column 95 Computing column 96 Computing column 97 Computing column 98 Computing column 99 Computing column 100 Computing column 101 Computing column 102 Computing column 103 Computing column 104 Computing column 105 Computing column 106 Computing column 107 Computing column 108 Computing column 109 Computing column 110 Computing column 111 Computing column 112 Computing column 113 Computing column 114 Computing column 115 Computing column 116 Computing column 117 Computing column 118 Computing column 119 Computing column 120 Plotting